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A continuum theory of partially fluidized granular flows is developed. The theory is based 
. on a combination of the equations for the flow velocity and shear stresses coupled with the order 

parameter equation which describes the transition between flowing and static components of the 
granular system. We apply this theory to several important granular problems: avalanche flow in 
' deep and shallow inclined layers, rotating drums and shear granular flows between two plates. We 

carry out quantitative comparisons between the theory and experiment. 

on ; 

PACS: 47.54.+r, 47.35.+i,46.10.+z,83.70.Fn 

o 

<D ■ I. INTRODUCTION 

The dynamics of granular material under shear stresses plays a fundamental role in many natural phenomena and 
technological applications. [Q-Q. When the shear stress exceeds certain threshold, granular material undergoes a 
transition from a solid state to a fluidized state (yield). The physical mechanism and properties of this transition are 
+J , still not completely understood. In many important situations, the granular material remains in a multi-phase state, 
when part of it is fluidized while another part is solid. 

On the theoretical side, a significant progress had been achieved by large-scale molecular dynamics simulations 
1 , ■ ^,[1 and by continuum theory |7|-p^| . The current continuum theory of dense near-surface flows was pioneered by 
Bouchaud, Cates, Ravi Prakash and Edwards (BCRE) || and subsequently developed by de Gennes, Boutreux and 
and Raphael 0Jl^,|ll[. In their model, the granular system is spatially separated into two phases, static and rolling. 
The interaction between the phases is implemented through certain conversion rates. This model describes certain 
features of thin near-surface granular flows including avalanches. However, due to its intrinsic assumptions, it only 
works when the granular material is well separated in a thin surface flow and an immobile bulk. In many practically 
, important situations, this distinction between "liquid" and "solid" phases is more subtle and itself is controlled by 
00 ' the dynamics. 

\ The purpose of this paper is to develop a unifying description of such partially fluidized granular flows and apply 
this theory to several problems of granular dynamics fl3|| . The underlying idea of our approach is borrowed from the 
Landau theory of phase transitions |Q. We assume that the shear stresses in a partially fluidized granular matter 
are composed of two parts: the dynamic part proportional to the shear strain rate, and the strain-independent (or 
"static") part. The relative magnitude of the static shear stress is controlled by the order parameter (OP) which 
varies from in the "liquid" phase to 1 in the "solid" phase. A possibility of describing a granular flow as a multi- 
phase system undergoing a phase transition, has been proposed by de Gennes without further elaboration. Unlike 
ordinary matter, the phase transition in granular matter is controlled not by the temperature, but the dynamics 
stresses themselves. In particular, the Mohr-Coloumb yield failure condition Q is equivalent to a critical melting 
temperature of a solid. The OP can be related to the local entropy (and possibly density) JlfJ of the granular 
material. OP dynamics is then coupled to the hydrodynamic equation for the granular flow. 

We apply this theory to several cases of granular flows of considerable interest. First, we will focus on gravity driven 
free-surface granular flows which typically occur in shallow chutes, sandpiles, and rotating drums. The most famous 
form of such flows is an avalanche, and our theory yields a rather detailed description of the avalanche dynamics. 
Then we apply our theory to granular Couette flows induced in the bulk by a moving boundary. Our model captures 
important phenomenology observed in these systems experimentally ]l6|-p7|]. 

The structure of the article is the following. In Sec. II we describe general formulation of the partially fluidized 
granular flows. In Sec. Ill we focus on free surface flow problem on incline plane (chute flow). In this section 
we consider stability properties of stationary solutions, avalanches in shallow chutes, transitions from triangular to 
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uphill avalanches, and comparison with experimental results. In Sec. IV we study flow in deep layers. We illustrate 
application of this theory for two-dimensional rotating drum. We show that our model exhibit avalanche flow at low 
rotation rates and transition to steady flow for higher rotation rates. In Sec. V we extend our approach for shear 
granular flow and discuss connection with dry friction phenomena in granular systems. In Sec. VI we discuss various 
implications of our results. 
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II. GENERAL FORMULATION 



We base the continuum description of granular flows on the momentum conservation equation 

Dvi da q 
~Dt ~ ~dx 



Po-R7 = +Po.g s :, J = 1,2,3. (1) 
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where Vi are the components of velocity, p — const is the density of material (we set po = 1), g is acceleration 
of gravity, and D/Dt — d t + Vid Xi denotes material derivative, cr^ denotes components of stress tensor. Since the 
relative density fluctuations are small, we assume that the velocity obeys the incompressibility condition V • v = 0. 
Here we emphasize that changes in density in fact are very important, since the onset of flow in granular materials 
itself is related to dilatancy, or small decrease in density. However, the relative changes of density for dense flows 
are typically very small (below few percent) , and therefore the compressibility effects for hydrodynamics of granular 
flows are negligibly small. Still, these small variations in density can substantially affect transport coefficients and 
constitutive relations. In fact, our order parameter equation describes sensitive response of granular media to small 
changes in local ordering, and, consequently, small changes in density. 

Momentum conservation equation (|l|) has to be augmented by the appropriate boundary conditions (BC) . As usual, 
on solid walls we require no-slip conditions Vi — 0, and on free surfaces, the kinematic boundary condition is assumed, 

= Vn (2) 

Dt ' 

where £ is the displacement of the free surface, and v n is the normal velocity to the surface. 

The main difficulties in describing granular flows center around the constitutive relationships for stresses <7ij . These 
relationships differ drastically for flowing and static configurations of granular matter. For static regimes, the shear 
stresses are determined by the applied forces, whereas in flows the shear stresses are proportional to shear strain rates. 
The transition from one regime to another is controlled by so-called yield criteria, among which the most popular is 
Mohr-Coloumb criterion relating shear and normal stresses. The goal of our paper is to unify the description of these 
different regimes of granular dynamics within a single theory. The central conjecture of our theory is that in partially 
fluidizcd flows, some of the grains are involved in plastic motion, while others maintain prolonged static contacts 
with their neighbors. Accordingly, we write the stress tensor as a sum of the hydrodynamic part proportional to the 
flow strain rate eij, and the strain-independent part, afj, i.e. ct^ = e<j + erf^. We assume that diagonal elements 
of the tensor coincide with the corresponding components of the "true" static stress tensor of- for the immobile 
grain configuration in the same geometry, and the shear stresses are reduced by the value of the order parameter p 
characterizing the "phase state" of granular matter. Thus, we write the stress tensor in the form 

^ = "(^ + fe) + (p+(1 -' ,)a ^- (3) 

Here rj is the viscosity coefficient, 5ij is Kronccker symbol. In a static state, p = 1, o~ij — of-, Vi = 0, whereas in a 
fully fluidized state p — 0, and the shear stresses are simply proportional to the strain rates as in ordinary fluids. 

To complete the set of governing equations, we need to introduce constitutive relations for components of the static 
shear stress tensor <xy , as well as an equation for the order parameter p. The issue of constitutive relations in static 
granular configurations is rather complex and not completely understood |§,^8). It appears that in many cases, the 
constitutive relations are determined by the construction history p9|. Recent studies elucidated the fundamental role 
of the force chains networks in formation of the shear stress tensor pQ] . We will assume that for any given problem, 
the corresponding static constitutive relations have been specified. 

Since in dense granular flows the energy is rapidly dissipated due to inelastic collisions, we apply pure dissipative 
dynamics for the order parameter p, which can be derived from the "free-energy" type functional T: 



Dp _ 8T 
Hi ~ ~~5p' 



(4) 



We adopt the standard Landau form for T ~ J dr(D\Wp\ 2 + f(p,<fi)), which includes a "local potential energy" 
and the diffusive spatial coupling. The potential energy density f(p, 4>) should have extrema at p — and p = 1 
corresponding to uniform solid and liquid phases. According to the Mohr-Coulomb yield criterion for non-cohesive 
grains j| or its generalization [^o|, the static equilibrium failure and transition to flow is controlled by the value of 
the non-dimensional ratio (f> = max|<7j{ m /<7jj n |, where the maximum is sought over all possible orthogonal directions 
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n and m in the bulk of the granular material. We simply use this ratio as a parameter in the potential energy for the 
OP p. Without loss of generality, we write the equation for p: 



T^ = l 2 V 2 p-p(l~ P )F(p,<j ) ) (5) 

Here r and I are characteristic time and length correspondingly. One can expect that the length I is of the order of 
the grain size and the time r ~ r g , where r g = \JTJg is typical time between collisions in dense granular flow. Further, 
according to observations, there are two angles which characterize the fluidization transition in the bulk granular 
material, an internal friction angle tan -1 <j)\ such that if <f> < <f>i the static equilibrium is unstable, and the "dynamic 
repose angle" tan -1 <f> such that at (f> < 4>o, the "dynamic" phase p = 0, is unstable. Values of 0o and </>i depend 
on microscopic properties of the granular material, and they do not coincide in general. Typically there is a range in 
which both static and dynamics phases co-exist (this is related to the so-called Bagnold hysteresis @ ) . The simplest 
form of F(p, (jj) which satisfies these constraints, is F{p, (j>) = —p + S, where 

« = (^-$)/(^-$) (6) 
Here we use a square of <j> to avoid non-analytical behavior at a® z — 0. Rescaling t — ► t/r and X{ — ► Xi/l leads to 

^=V 2 P + P(1-P)(P~5). (7) 

This equation completes the general formulation of the continuum theory for partially fluidized granular flows. In 
an infinite system with fixed stress parameter < 6 < 1 (4>o < 4> < 4>i), both static (p — 1) and dynamic (p = 0) 
phases are linearly stable, and Eq. (|7j) possesses a moving front solution which "connects" these phases. The speed 
of the front in the direction of p = is given by V = (1 — 2S)/y/2. At 8 = 1/2 both phases co-exist. For 5 < 0, 
only "solid" phase survives, and at 5 > 1, only liquid phase survives. The dynamics of partially fluidized granular 
flows becomes much more interesting in confined systems with fixed or free boundaries. In the following sections we 
consider several such problems of particular interest. 



III. SHALLOW GRANULAR FLOW ON AN INCLINED PLANE 



Let us now specialize our theory to the description of free-surface chute flows. We consider a layer of dry cohesionless 
grains on a sticky surface tilted by angle if to the horizon. We introduce a Cartesian coordinate frame aligned with 
the tilted layer (see Figure Q). In case of a stationary shear flow, the force balance of Eq.([l]) yields the following 
conditions: 

&zz,z + <7xz,x = ~g COS If , (T xz ,z+ &xx,x = 9 SU1 if (8) 

where the subscripts after commas mean partial derivatives. The solution to Eqs. (||) in the absence of lateral stresses 
°yy = a yx = <?yz = 0, is given by 

ct Z2 = -gcospz , a xz = gsinipz , a xx , x = (9) 




FIG. 1. Schematic representation of a chute geometry 
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Thus, in a stationary flow there is a simple relation between shear and normal stresses, a xz — — tan tpo zz independent 
on the flow profile. In a static equilibrium, the force balance also gives a xz = — tan ycr° z . Since by assumption 
a zz = a zz , we obtain a xz = a xz . In a flowing regime, the total stresses are composed of the static contribution in 
the viscous strain-related terms. According to our conjecture, the same relation holds between the static parts of the 
stress in the flowing regime. In this Section, we will consider non-stationary process of avalanche flow, but we will 
assume that this simple constitutive relation between shear and normal stresses is maintained in this regime as well, 
and deviations from the stationary stress distribution are small. For the chute flow geometry, the value of parameter 
4> in Eq. (jfy can be easily specified. In this case, the most "unstable" yield direction is parallel to the inclined plane, 
so we can simply write <fi = \cr xz /a zz \ = tany>. On the free surface we impose the no-flux condition for the order 
parameter, p z = 0, and at the bottom z — —h we set p = 1 (a granular medium is in a solid phase near the no-slip 
surface). All components of velocity should be zero at the bottom z — —h. The kinematic boundary condition (^|) on 
the free surface for incompressible medium can be expressed in a form of the mass conservation law 

d t h = -(d x J x + d v J y ), (10) 

where J x ^ y — f , v x ^ y dz are in-plane components of the flux of the granular material. In a typical situation of the 
chute flows, the downhill velocity v x is much larger than the orthogonal y-component v y , so the mass conservation 
constraint can be simply expressed as 

d t h = —d x J. (11) 

The velocity v x is determined from the order parameter via Eq.(|3|) with the no-slip boundary condition v x = at 
z = —h. 

The mass conservation law Eq. (^) can be rewritten in term of the parameter 5 which is related to the local slope 
d x h — <j). If we assume that the difference between critical values 01,2 is small, i.e. (<pi — 4>q)/4>i <C 1, which is the 
case for most granular flows, and \d x h\ <C tany, from Eq. @ we obtain (see also Fig. Q) 

<t> = d x htt~(6-6 ) (12) 

where (3 = 1/ {4>i — 4>q) > and S = const corresponds to unperturbed value of 8, i.e. describes the flow with constant 
thickness h. Substituting Eq. ( |l2| ) into Eq. ([ll]) one derives 

d t S = (3d 2 x J (13) 

This equation should be used instead of the conservation law ([ll]) for infinitely deep layers (sandpiles or heaps), where 
the thickness h is not defined. 

Because of the no-slip boundary condition at the bottom of the chute , for s hallow layers the flow velocity is small, 



so the convective flux of the order parameter can be neglected (see Sec. [II B for details), and the material derivative 
Dp/Dt in Eq.(0) can be replaced by dtp, 

d tP = V 2 p + p{l- p)(p-5). (14) 



However, for fast flows this term may become important, see discussion below in Sec. [VB 



A. Stationary solutions and their stability 

There always exists a stationary solution to Eq. ( |l4| ) p = 1 corresponding to a static equilibrium. For 5 > 1 it 
is stable at small h, but loses stability at a certain threshold h c > 1. The most "dangerous" mode of instability 
satisfying the above boundary conditions, is of the form p — 1 — Ae xt cos(7rz/2/i), A <C 1. The eigenvalue of this mode 
is 

X(h) =6-1- tt 2 /Ah 2 (15) 
Hence the neutral curve A = for the linear stability of the solution p = 1 is given by 

For h > h c {8) grains spontaneously start to roll, and a granular flow ensues. 
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In addition to the trivial state p = 1, for large enough h there exists non-trivial stationary solutions satisfying 
the above BC. These solutions correspond to bold lines in the phase plane plots Fig. 0. These solutions describe 
stationary granular flows supported by a constant supply of granular material up-stream. For 1/2 < S < 1 there exists 
a separatrix of the saddle p — 1, p z — corresponds to the localized near-surface flow in infinitely deep layer. 



8<0.5 




1\ P 




FIG. 2. Phase plane of stationary Eq.(|14[) for three typical values of 8: a, 5 < 0.5; b, 0.5 < S < 1; c, <5 > 1 

The velocity profile corresponding to a stationary profile of p(z), can be easily found from Eq. (||) taking into 
account that a xz = a®, 



-^j = (1 - p)(T^ = -fl(l - p)z, 

where p = g sin <p/r). The flux of grains in the stationary flow J is given by 

v x {z)dz = —p / / (1 — p(z'))z' dz' dz = p / z 2 (l — p)rfz. 



(17) 



(18) 
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FIG. 3. Thickness of the stationary flow h as a function of supplied granular flux J/p, at several 5 (the value of 6 is shown 
near each line). The parts of the curves corresponding to unstable solutions, are shown by dashed lines. The thick line, given 
by the envelope condition dJ(h, 5)/dS = 0, separates stable and unstable parts of the upper branches. Thick dashed line shows 
the stability exchange line dJ(h, 8)/dh — 0. Inset: Normalized slope S vs flux J/p for infinitely deep chute (h — > oo), the OP 
p is given by Eq. (pi]). The branch below the dashed line (the condition dj(6)/d8 = 0) is unstable. 

The flux of supplied granular material J controls the thickness of the layer and the velocity profile. Figure || shows 
the thickness of the layer as a function of flux J at several values of S. For a fixed J, there are two values of h 
which correspond to two different regimes of granular flow on an inclined plane. The lower branch corresponds to 
the flow which involves the whole layer, while the upper branch corresponds to a flow which is localized near the 
surface. The selection and the stability of these solutions depend strongly on the particular problem at hand. Since 
for fixed h Eq. ((ij) has a free energy functional (see Eq (£|), T = \(p\ + p 4 /2 + 5p 2 - 2/3(6 + l)p 3 ), stable solutions 
should correspond to the minimum of T . It is easy to check that the lower branch corresponds to minima of the free 
energy T, and therefore it is stable, whereas the upper branch is unstable (corresponds to the maximum of J-). If 
the flux J is fixed by the boundary condition at x = 0, the corresponding unstable mode, which at large h is close 
to the translation mode d z p, is prohibited because it would violate the mass conservation constraint. However, even 
in this case the solutions corresponding to the upper branch can be unstable with respect to spatially non-uniform 
perturbations. Actually, it can be easily shown that this type of instability occurs for parts of the upper branch 
solutions corresponding to dJ/dS < 0. Indeed, assuming small deviations of the local slope Si with respect to average 
5 and linearizing Eq. (pa), we obtain 



9tSi 



(19) 



For dJ/d5 < it is a diffusion equation with negative diffusion coefficient, which is subject to a long-wavelength 
instability. On the other hand, if the slope of the free surface is fixed everywhere, then this instability is also 
suppressed. 

The two branches merge at the minimum value h s (S) where dJ/dh = 0. At h < h s , there is no stationary granular 
flow solution, and only non-stationary regimes are possible (see below). The value of h s can be found as a minimum 
of the following integral as a function of po, the value of p at the surface z = 0, 



h. = 



dp 



P0 ^.mpL + sf-cfa) 



(20) 



where c(po) — Pq/2 — 2(S+1)pI/3 + Spq. This integral can be calculated analytically for 6 — > oo and 6 — > 1/2. It is easy 
to show that for large 6, the critical solution of Eq.(|l4|) has a form p = 1 + Acos(kz) with Ac 1 and k = (6 — l) 1 / 2 , 
and therefore, h s (5) — > h c (5). For 5 — > 1/2, the critical phase trajectory comes close to two saddle points p = and 
p = 1, and an asymptotic evaluation of ( p0| ) gives h s = -V21og(£- 1/21+ const. This expression agrees qualitatively 
with the empirical formula <f> — (j>o ~ exp[—h s /ho] proposed in Ref. |19|j2l[| . 
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Neutral stability curve h c {5) and the critical line h s {8) limiting the region of existence of non-trivial stationary 
granular flow solutions, are shown in Fig. Q They divide the parameter plane (8, h) in three regions. At ft, < h s (S), 
the trivial static equilibrium p — 1 is the only stationary solution of Eq.(|l4|) for chosen BC. For h s (6) < h < h c (8), 
there is a bistable regime, the static equilibrium state co-exists with the stationary flow. For h > h c (S), the static 
regime is linearly unstable, and the only stable regime corresponds to the granular flow. This qualitative picture 
completely agrees with the recent experimental findings |ll],|2l| . Moreover, for no-slip bottom BC (corresponding to 
our p=l), authors of Ref. found a region of bistability in the parameter plane {h,(p) which has a shape very 
similar to our stability diagram Fig. | (see below Sec Jm^). 



25 



20 



15 



10 



solid & 
flow 



solid only 





1.08 1.13 1.18 1.23 

5 

B flow only 

quasi-linear limit 



0.25 



0.75 



1.25 



1.75 



FIG. 4. Stability diagram. Dashed line shows neutral curve (0), solid line shows the existence limit of fluidized state (po|). 
Dot-dashed curve depicts the transition line from triangular to uphill avalanches obtained from solution of Eqs. (|l4|),(|ll"|) for 
fi = 0.4 and (3 = 0.25. The line with circles shows the results obtained in quasi-linear limit, Eqs. (pi|),(fj5|) for j3 = 0.25 and 
the same value of [i which corresponds toa = 0.05. Inset: ip (in degrees) vs 8 for f3 = 0.25 and a — 0.015. 



As we mentioned above, the upper branches of the h(J) curves in Fig.^ correspond to the case of a near-surface 
flow. For large enough h this regime can become unstable with respect to spontaneous change of the slope 8. As was 
outlined above, the change of stability occurs at a tangent point between a curve h(J) and an envelope h e (J) to the 
family of curves h( J) for various 6, where 9,5 J = 0. The instability would exhibit itself as accumulation of granular 
material near the top of the inclined plane leading to the change of slope. This process will result in an unlimited 
growth of local depth h, and at t — > oo, the new stationary solution corresponding to h — > oo will be achieved. This 
regime can be described by an analytical formula which corresponds to the separatrix in Fig. ||b, (cf. Ref. |3l[| ): 



P = 



y/(S + l){8- l/2)cosh(zVT _ 8) + 28-1 
y/(S + 1) (8 - 1/2) cosh(z VT~S) + 2-8 



(21) 



In this deep-layer solution, the parameter 8 which corresponds to the slope of the free surface, is not related to the 
slope of the inclined plane (free surface can be more or less steep than the underlying plane, as in sandpiles). Rather, 
8 is determined by the value of J. The dependence of the slope 8 vs. flux J for solution (|2l]) is shown in Fig. ||,inset. 
The condition J = const gives rise to two stationary values of 8. The upper branch approaches 8 — 1 as J — * oo as 
8. For the lower branch, the width of the fluidized zone zq, defined by p(z = zq) = 1/2 is growing as 



J r 



1/VT 
log(S 



1/2) for 8 — ► 1/2. Correspondingly, in this case one has the relation between the flux J vs 8: 



J 



z 2 {l- p)dz~ zl~ |log(5-l/2)| 3 



(22) 



Both branches merge at some minimum J = J c . In the vicinity of J c the flux and the angle are related as J ~ 
J c + const x (8 — 8 C ) 2 + .... According to the condition dJ/88 > 0, only the upper one corresponds to a stable 
near-surface flow, and the lower one corresponds to an unstable regime. In the stable regime, the slope of the sandpilc 
increases with the flow, and for very large J the slope of the free surface 8 approaches 1. This behavior agrees 
qualitatively with observations Ref. [£7j . 
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However, if the change of 8 is constrained, the instability of the lower branch can be suppressed. We believe that 
in Ref. [^6| , placing the mouth of the hopper supplying the sand directly on the surface of the sandpile limits the 
variation of 6 and may possibly stabilize the lower branch. Moreover, since the instability is of convective type, the 
length of the system may not be sufficient for it to develop. 

At J < J c , a stationary flow does not exist. In this regime, the granular material accumulates and discharges in a 
form of avalanches periodically in time (see below). This phenomenology is also consistent with recent experiments 
in Ref. [ p7|]35| ] where the transition from intermediate avalanches to steady flow is reported. Moreover, as one 
expects from the mass conservation, if the flow is represented by periodic sequence of well-separated avalanches each 
carrying an amount of grain, the time between the consequent avalanches should be Tq ~ J \ in agreement with 
experiment p7[. O ur theory predicts that at the onset of steady flow the angle of sandpile should show critical behavior 
8 — 8 C ~ y 'J — J c . However, in experiment p7| the critical transition has not been detected, possibly because the 
slope changes within a narrow range. 

Ref. j27| also reports an increase of the width of fluidized layer zq with increase of the applied flux J, which is 
consistent with Eq. Also in agreement with the theory, the heap angle in Ref. |27j increased with J, which 

corresponds to the upper (stable) branch of the dependence 8(J) shown in Fig. |],inset. 

Similar transition is known for partially-filled rotating drums when the rotation speed is varied. For low rotation 
speeds the flow in the drum occurs in the form of periodic sequence of avalan ches, whereas for larger rotation speeds 
a steady surface flow ensues |52| . We discuss this case below in Section [V C . 

We compared the velocity profiles measured in Ref s. [f27| , p6[ with our theory. The velocity can be determined from 
Eq. ([IT]) using expression for the order parameter (|2l|). A typical velocity profile v(z) vs z is shown in Fig. 0. 
For convenience we scaled v(z) by the value at open surface v(0). In agreement with the experimental data, the 
stationary profile has an exponential tail, i.e v(z) ~ exp(— z/d s ) where d s = — 8. For the lower branch of S(J) 

which apparently describes the flow in this experiment by Komatsu et al. Eq), the parameter 5 at high flow rate 
approaches 1/2 (see inset of Fig.||), i.e. the decay length d s is d s = 1^/2 w Z/0.707. Here I is the characteristic length 
in Eq. (||). Experimental value is d s w d/0.72, which fixes the characteristic length equal to the grain size, i.e. I = d. 
Moreover, experimental data of Ref. |26| strongly indicates independence of the decay length d s on the value of flux 
in wide range of flux values and grain diameter. This behavior is again corresponds to the lower branch of 8{J) 
dependence. The value of the characteristic length / agrees with other independent experimental observations (see, 
for example, Refs. f p^pil) . 

Lameux and Durian |27|] also found exponential decay of velocity down the surface: v ~ exp[— z/(0.15cm)] for the 
grain diameter d = 0.33 ± 0.03 mm. In their experiment, unlike Ref. [^6), the particles were allowed to fall on the 
top of the sandpile, thereby relaxing the constraint on the slop of the sandpile. In this case, an upper branch of the 
function 8 (J) should be selected, and for that branch, the slope 8 > 0.72, so the characteristic decay length indeed 
should be larger than particle size d. In fact, it should be directly proportional to the flux J. It would be interesting 
to test this prediction in future experiments. 




FIG. 5. Stationary velocity profiles viz) vs the distance from free surface z/d for different grain sizes d (mm) and supplied 
flux values J (grams/sec), from Ref. M. Solid lines show theoretical results for two different values of S. 
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Pouliquen [Ell proposed a scaling for the mean velocity v — J/h vs thickness of the layer h in the stationary flow 
regime, v oc hr^/h s , which works for various angles ip as well as for different grain sizes. Eq. (|l§| ) yields v oc (h — hs) 1 ^ 2 
for small h — h s and v oc h 2 for large h. It is plausible that the experimentally found scaling exponent 3/2 is the 
result of the crossover between the two different regimes. However, renormalization v/\fgh, h/h s as in Ref. H| does 
not collapse our results onto a single curve, perhaps due to the assumption of a simple Newtonian relation between 
the strain v z and the hydrodynamic part of the shear stress a xz with a fixed viscosity rj (see Eq.(||)). In fact, the 
viscosity itself may depend on p and z in some fashion. 



B. Non-stationary dynamics in a single mode approximation 



form (compare Sec. |III A| ) 



In the vicinity of the n eutral curve (fLq) Eqs.(0lJl4) can be significantly simplified. We may look for solution in the 



p= 1- A(x,y,t) cob (j^z) +0(A 2 ), (23) 

where A <C 1 is now a slowly varying function of t, x, and y. At the neutral curve defined by the condition 
X{5, h) = 5—1 — tt 2 /Ah 2 = the expression ( p3| ) with A = const, h = const is an exact solution to linearized Eq. (fil|). 
In the vicinity of the neutral curve defined by condition |A| <C 1 the ansatz ( p3| ) with the slowly- varying functions 
A, h gives an approximate solution to full Eq. ([l4]). The function A itself is determined as a result of orthogonality 
(or solvability) condition with respect to function cos (gr^)- 

Substituting ansatz ( p3|) into Eq. (Q) and applying orthogonality conditions, we obtain in the first order 

A t = XA + V\A + W-^ A 2 - - A 3 - ah 2 Ad x A (24) 

where = d 2 + d 2 , and a = p(3ir 2 — 16)/37r 3 = 0.146/j. Eq. (]24|) must be coupled to the mass conservation 
equations which reads as (here we neglect contribution from the flux along y-axis J y ~ d y h <§; J): 

dh _ dJ _ a 9h 3 A 
dt dx dx ' 

where J was calculated from Eq. (pit) using ansatz (p3|) and a = 2fi(ir 2 — 8)/ir 3 = 0.12/i. Taking into account that 



variations in h also change local surface slope, we replace 5 in (24) by 5q — f3h x , see Eq. (|12|) 

In deriving this equations we assumed that (2 — 5) A 2 and A 3 are of the same order, i.e. S w 2, however qualitatively 
similar equation with a different nonlinearity can be obtained for any S and h. 

The last term in Eq. (^) originates from the convective term vV p in Eq. ([?]). For not very large thicknesses of the 
layer h and in the large viscosity limit (which we actually consider in the paper) p <^ 1 this term can be neglected 
with respect to other terms in Eq. (p4j). However, for thick layers the convective term cannot be neglected because 
the magnitude of this term grows as hr. 

We studied Eqs. (p4|) , (p5|) numerically. First, we considered the flow with the fixed supplied flux in one dimension 
(x). In this situation the flux J = const is introduced via the boundary condition in Eq. ( |25| ) at x — 0. For large 
values of the flux we indeed observed the transition to the steady flux regime (see Fig. ^a), although some transient 
avalanches occur which are related to the adjustment of the chute thickness. For smaller values of the flux (below the 
corresponding cutoff value on Fig. |^) , we find that the flow occurs in the form of periodic sequence of avalanches, Fig. 
^b. Our numerical simulations indicate that the time between the avalanches To diverges as Tq ~ J -1 at J — > 0, in 
agreement with experimental result in Ref. J^?]]. Moreover, we observed abrupt hysteretic transition from avalanching 
to steady flow with the increase of supplied flux J, which also agrees with Ref. p7| . 

In order to study the evolution of avalanches in two dimensions (a;, y) we performed simulations in a fairly large 
system, L x = 400 dimensionless units in x-direction (downhill), and L y = 200 units in y-direction, with the number of 
grid points 1200 x 600 correspondingly. As initial conditions we used uniform static layer: h = ho, A = 0. We triggered 
avalanches by a localized perturbation introduced near the point (x, y) = [L x /A, L y /2). Close to the solid line in Fig. 
(0) we indeed observed avalanches propagating only downhill, with the shape very similar to the experimental one 
p0| . The avalanche leaves triangular track with the opening angle ip in which the layer thickness h is decreased with 
respect to the original value h$. At the front of the avalanche the layer depth is greater than ho, as in experiment. 
The opening angle as a function of 8 is shown in inset of Fig. 0. 
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FIG. 6. Space-time surfaces showing the ID evolution of height h in a shallow chute for two values of the fixed supplied 
flux J = 0.6 (a) and J = 0.4 (b). Other parameters of the model are a — 0.025, /3 = 3.15,5 = 1. The chute length L = 500, 
number of grid points N = 1000. Initial condition is A = 0, h — 3. 




FIG. 7. Grey-coded images demonstrating the evolution of a triangular avalanche for t = 50 (a), t = 200 (b) and 250 (c). 
White shade correspond to maximum height of the layer, and black to minimum height. Parameters of Eqs. (pi|) , (pBj) are: 
a = 0.15, P = 0.25, S = 1.2 and h = 3, point A in Fig. g. 

For larger values of 5 or for thicker layers (close to the dashed line in Fig. ||) we observed avalanches of the second 
type. In this case the avalanche propagates also uphill, and, unlike the previous case, the whole avalanche zone is in 
motion, as new rolling particles are constantly arrive from the upper boundary of the avalanche zone. Sometimes we 
observed small secondary avalanches in the wake of large primary avalanche, see Fig. ||c. 



10 



FIG. 8. Snapshots of an uphill 




arc: 



a = 0.05, /3 = 0.25, 8 = 1.07 and ho = 5.5, point B in Fig. W. A small secondary avalanche is seen on the image (c). 



C. Transition from triangular to uphill avalanches 



Our model predicts the transition from triangular to uphill avalanches when the thickness of the layer or the 
inclination angle are increased, similar to that observed in experiment Jl9[ |. In order to investigate the transition in 
detail, we again return to the one-dimensional version of Eqs. (pq). The result of simulations are shown in Fig. 
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FIG. 9. ID evolution of a localized perturbation in a long shallow chute for two values of S, (a) 5 = 1.02; (b) S — 1.07. 
Shown are the height profiles at 10 consecutive moments of time for ho = 5.5, a — 0.05, (3 = 0.25. Secondary avalanches are 
seen on panel (a). 

As can be seen in Figure ^]a, finite perturbation introduced near x = 300 triggers a downhill avalanche for smaller 
5. Due to mass conservation the height of the avalanche increases as it propagates downhill. 

For larger S, the region of fluidized grains grows not only downhill, but also uphill (see Fig. [|b). In contrast to 
the downhill avalanche, the uphill front appears to be a steady-state solution, A = A(x + Vt), h = h(x + Vt). Our 
simulations show that the velocity of the front remains finite at the transition point, see Fig. mi. Since the uphill 
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front always propagates with the velocity V > Vq > 0, we call this phenomenon "velocity gap". 




1.10 



FIG. 10. The velocity of uphill front vs 5 for ho = 5.5, a = 0.05, j3 = 0.25. The cutoff value of velocity Vo w 0.35. 



It can be shown from the analysis of Eq. ( pq ) that the front solution indeed cannot exist for arbitrarily small speed. 
In a co-moving with velocity V frame the front is a stationary solution, and Eq. (feq) reads 



V(h-ho) = 



*h 3 A 



The non-trivial front solution must satisfy the boundary conditions h 
Aco 7^ for x — > co, where 



ho, A — > for x 



-co and h 



(26) 
,A = 



A, 



9tt 




16(2-5) {16(2 -6) 



97T 



rU-i- 



4A§ 



(27) 



Since Aoc, cannot be arbitrarily small for finite hoc, by the nature of the hysteretic transition from solid to the fluidized 
state, V also cannot be arbitrarily small. Thus, our model predicts the velocity gap for the uphill front, which is in 
fact supported by the experimental data in Ref. [2(J . This result appears to be in contradiction with the conjecture 
of Bouchaud and Cates |$4j that the transition from triangular to uphill avalanches occurs at zero front velocity. 

Tracking systematically the moving front existence limit in (5, h) we obtained the line separating the triangu- 
lar/uphill avalanches in (6,h) plane, see Fig. |]. 



D. Uphill-triangular avalanche transition and the velocity gap in the large viscosity limit 

The velocity gap disappears in the infinite viscosity limit (i.e. a = 0). In this case, the thickness of the layer does 
not change (h = ho = const), and the order parameter equation ( |l4| ) becomes independent. The uphill front solution 
p(x + Vt, z) satisfies the equation 

V Px = V 2 p-p(l-p)(5-p) (28) 
In this case, the transition between uphill and downhill front propagation is continuous, and can be obtained from 



the solution to Eq.(28) with V = 0. This solution exists only for a specific value of 5 corresponding to the layer 
thickness h. The dependence 5(h) is derived in Appendix |A| At large ho it results in 

log(*-l/2) 

h 5-1/2 ' 1 yj 

i.e. at large h the region of uphill avalanches shrinks, in agreement with experiments fl9|| . 

For small h = 0(1), the transition line can be found from the stationary solution of single-mode approximation 
Eq. (p4|) . In this case, 
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h=\ , C'.O! 



(see Appendix |X|) . 

For small but finite a (large viscosity rj), the velocity gap is small (C^a 1 / 2 )), and can also be found analytically 
from Eqs. @, @. 

At a <C 1, the uphill front speed satisfies the following equation (see Appendix |b]) 

7_ Wl + ^ 2 -=0, (31) 

where 8 = 8 — 8* and 8* is determined from V = condition at a = 0, and di^ are specified in Appendix |^. From 
Eq. (fl|) we find 



V = ^Y + ^]{di~5) 2 /±-ad 2 (32) 

(the branch with "— " sign in front of the square root is unstable). Thus, the cutoff value of the velocity Vq = ^ and 
corresponding value of 8 at the threshold of uphill propagation is 

8 = 6* + 2 v ^h/d 1 (33) 
The above expansion however is valid only for very small a obeying the condition ah 3 <C 1. 

E. Comparison with experiment 

In order to establish a link between our theory and the experiments we need to specify the parameters 4>o and </>i , 
as well as characteristic length I and time r, and the viscosity rj. Parameter 4>i can be easily determined from the 
value of the chute angle corresponding to the vertical asymptote of the stability curve on the experimental bifurcation 
diagram of Ref. M]. The value of (f>o cannot be directly read from the bifurcation diagram. However, the vertical 
asymptote to the line bounding to the region of existence of avalanches in Ref. JlSj ] , gives the value of the angle 
(j)Q which corresponds to the regime of when front between granular solid and fluid does not move, i.e. 6=1/2 
instead of 6 = 0. Thus we can express our parameter 8 through cf>o,(j>i. For the experimental parameters of Ref. |fl9| , 
tan -1 0o ~ 25° and tan -1 <f>i w 32°. It gives (3 = 1/2(01 — 0o) ~ 3.15. Based on the comparison with experimental 
results for velocity decay in stationary flow Ref. Pq| , as a characteristic length I we can take the mean diameter of 
the grain d which for experiment Ref. jl9| was 0.24 mm. Solid and dashed lines in Fig.[ll] indicate theoretical stability 
boundaries, which correspond very nicely to the experimental findings. 

The position of the line separating the triangular and uphill avalanches depends on the value of parameter a in Eq. 
(^B|). In fact, a ~ r/rj is the only fitting parameter in the theory. In principle, it could be determined independently 
if we knew the characteristic time and the viscosity, but this data is not available to us. We find from numerical 



solution of Eq. (|24|), (25) that the best fit to experimental data occurs for a ~ 0.025 (correspondingly /i w 0.2). For 
this choice one observes a very good correspondence between theory and experiment (dotted line in Fig. |lT| ) . 
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FIG. 11. Comparison of theoretical and experimental phase diagrams. Lines obtained from theory, symbols depicts exper- 
imental data from Ref. Jl^]. Solid line and circles limit the range of existence of avalanches, long-dashed line and triangles 
correspond to the linear stability boundary of the static chute, and the dot-dashed line and crosses denote the boundary between 
triangular and uphill avalanches for f3 = 3.15, a = 0.025 (or, correspondingly, p — 0.2). Dotted line shows the transition curve 
in the infinite viscosity limit r\ — > oo. 



IV. FLOWS IN DEEP GRANULAR LAYERS 



A. Avalanches in deep chutes 

In deep granular layers, our assumption that the convective flux of the order parameter is small, is no longer 
valid, and we have to return to Eq.(^). For smooth horizontal variations of the flow, its local vertical profile can be 
approximated by the following dependence at — oo < z < 0, 

p = l- (tanh[(z + z )/VS] - tanh[(z - z Q ) /Vsfj /2. (34) 

with slowly varying depth of the fiuidized layer z$. This expression is very close to the exact front solution 

- (l ±tanh[z/V8]) (35) 



9 2 , 

if zq 3> 1 and 8 — ► 1/2 and differs from it only in the vicinity of the free surface z = where it is augmented in 
order to satisfy the no-flux boundary condition d z p = 0. Moreover, for zq — > one has p — > 1, thus one recovers the 
behavior of linearized Eq. (Q) . 
Let us introduce the new variable 

o 

(1 - p)dz. (36) 



It is easy to check that for ansatz (34) in fact z = Zq. We will show below that the simplified description of the 
dynamics of Eq. (]?]) in the framework of z is rigorous in two important limits: z> 1 and 2<1. For the intermediate 
values of z the above approximation for the order parameter Eq. ([34]) gives smooth interpolation between these 
two limits. Our numerical simulations indicate that qualitative features are not sensitive to the specific choice of 
interpolation since the solution tends to "avoid" the intermediate area (we obtained qualitatively similar results using 
piece- linear approximations). 

After integration of Eq. (^) we obtain 

d t z = d 2 x z+ / p(l - p)(5 - p)dz + / (v x d x p + v z d z p)dz. (37) 



Horizontal velocity profile v x (z) is found from Eq. (p"7j) , 

v x = -p. [ (1 - p)z'dz'. (38) 



dz'd x v x = d x z p dz' dCC^cl 1 ~~ P)- ( 39 ) 



and 



Now, substituting Eqs. (|34D,(|3q),(|39D into Eq. (|37j), after some algebra we get 

d t z = d x z + F(z ) - pG(z )d x z . (40) 
Function F (zrj) can be found in the closed form 

F= « - + ^-*«L(* +S+1 ) (41) 



V2(s - 1) V2 
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with s = exp(-\/2zo). Function F(zq) has the following asymptotic behaviors 



F(z ) 



(5 — l)z for z <C 1 
V2 (6 - i) for z > 1 



(42) 



Thus, at small zq Eq. ( |37| ) complies with the behavior of the linearized near p = 1 Eq. (|7j), and for large zo and Eq. 
([37]) gives the asymptotically correct result for the velocity of the front between fluidized and solid state at S — ► 1/2. 

Function G(zo) can only be found in an integral form. However, asymptotic values of G(zq) can be found for large 
and small zq, 



G(z ) 



12-7T z 

3V2 



z sa 0.5021z for z < 1 



3.29 



3 ~ for z > 1 

The expression for G, valid also for intermediate values of zq, can be approximated as 

- 2 '12-ir 2 



G(z a ) = y tanh 



tt 2 V2 



■-0 



(43) 



(44) 



This equation has to be solved together with the equation for <5. The latter can be derived from the the mass 
conservation Eq. (|ll]). with the expression for flux given by Eq. (|l8|). Substituting p(z) from (|34|), we obtain 



where 



/(*>) 



dh |J„ > 
^ = -3^0), 



27r 2 z for zq <C 1 



for z > 1 



(45) 



(46) 



We used the simplest interpolation for this function: f(zo) = zq(z 2 
we arrive at the equation for S (compare Eq. (Jl3])), 

dtS - ffd 2 J(z ), 



2tt 2 ). Differentiating Eq.(^5|) with respect to x, 



(47) 



Equations (|40|),(47) give a simplified description of two-dimensional flows in deep inclined layers or sandpiles. We 
performed numerical simulations of this model in application to avalanches. We have found that small localized 
perturbations decay, and large enough perturbations trigger an avalanche. Figure [l^ shows the development of the 
avalanche from a localized perturbation imposed at the point x = 480. As it is seen from the Figure, the avalanche 
propagates both uphill and downhill. This observation is consistent with our conclusion from previous sections that 
the domain of existence of triangular avalanches shrinks with the increase of layer thickness. 



initial perturbation 




1 ' ' • 1 ■ ' ' 1 

200 400 600 800 

X. 

FIG. 12. Evolution of a free surface profile during an avalanche within a simplified model (^o|),(^7|) for 

r5n = 0.75, fi — 0.2, (3 — 3.15, i.e. the parameters are the same as for Fig. [H]. In the wake of the avalanche the slope of 
the free surface is reduced and approaches the equilibrium value 1/2. Note the "true" horizontal and vertical variables (a;,, z„) 
which are related to our original Cartesian variables (x,y) via a simple rotation transformation by angle ip. 
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B. Connection with BCRE theory 



It is interesting to point out the similarities and differences between our Eqs. ( 40 ) , (ft5|) and the set of phenomeno- 
logical equations for avalanches in deep layers proposed earlier in Refs. |f7|-|l0||. 

The BCRE theory || operates with two variables, thickness of immobile fraction H and the thickness of rolling 
(flowing) fraction R. These quantities obey the set of equations (see, e.g. Ref. 0) 

9R fl dR d 2 R 

- = -^ r -VR-,- + D — (48) 

— = 7 (0r - <P)R (49) 

where </> r is critical slope, 4> = —dH/dx is the local slope, term j((f>r — 4>)R describes the mass exchange between 
rolling (R) and static (H) layers, v is the flow velocity (assumed to be constant) within the rolling layer, and D is the 
diffusion constant. The first BCRE equation (4S) describes the dynamics of the rolling fraction, and equation ( |49| ) 
is analogous to our mass conservation law Eq.(45), although in our description h indicates the total thickness of the 
layer, i.e. h = H + R. 

These equations were later modified in Ref. []lo|-|l2f for flows involving large values of R by replacing the "instability 
term" 7(0 — (j) r )R by the "saturation term" (<f> — (j) r )v up for R > i? , Ro ^ 1, yielding 

OR . 1 \ dR d 2 R 

^ = ^~^)v up + v- + D^ (50) 

where v up is the constant of the order of v. This modification provides layer thickness saturation at large R. 

One may notice that Eq. ( |40| ) with ([H]),(Q) coincide with the first BCRE equation for zo <C 1, and, respectively, 
Eq. ( p0| ) for large Zq 7 however with one important caveat. From our derivation it directly follows that the value of the 
critical angle cf> r must be different for small and large R, whereas in Eqs.(f48|),(|5"c|) that value is kept the same. This 
important distinction of our model gives rise to the hysteretic behavior of the fluidization transition which is missing 
in the original BCRE model and its later modification. 

In their recent work Aradian, Raphael and de Gennes |l2| added phenomenologically the dependence of the velocity 
profile on the flowing layer thickness R. Note that in our approach this dependence appears naturally, however the 



particular form of the the coefficient G at the convective term in the equation (40) differs from a simple linear form 
proposed in Ref. |12|. 



C. Flow in rotating drum 



The dynamics of granular material placed in a rotating drum is another example of partially fluidized granular flow 
in a deep granular layer, see for review ||32| . Depending on the rotation rate u), see Fig. hj, the flow occurs in the 
form of sequence of avalanches for small u>, or steady flow for larger rotation speeds. At relatively small ui, the flow is 
confined to a narrow near-surface region, and the bulk exhibits rigid body rotation. These observations prompted a 
number of recently introduced continuum models ||7|, |l0| , ^8U4l| , ^1 in which the flow was described by the dependence 
of its total flux on the local free surface slope S, J ~ tan<5 — 8 C , where 5 C is the critical slope. This description 
yields rather realistic profiles of the free surface in the stationary flow regime, and also can explain main features 
of segregation of binary granular mixtures in rotating drums [|38|-|42]. However, it fails to describe the transition 
from periodic avalanching to the stationary regime. We believe that this can only be done within a model which 
incorporates the hysteretic character of granular fluidization. 



1G 
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FIG. 13. Sketch of a flow in a rotating drum 

In this section we focus on the non-stationary granular flows in 2D rotating drums, see Fig. [l^. To simplify the 
description, we assume that the free surface profile is not very different from a straight line and the radius of the 



drum is much larger then the grain size. It allows us to use the formalism developed in Sec. IV A. We reduce the 
description of the flow to evolution of only two quantities: the position of the solid/liquid interface zq and parameter 
6 which is proportional to tan 0, i.e. local slope of the free surface. 



Equation for zq is the same as in Sec. IV A 



d t z = d s z a + F(z ,S) - vd s z , (51) 

where we introduced the the coordinate s along the free surface and the convective term vd s Zo, see discussion above. 
This equation is subject to boundary conditions zq = at the drum walls (say at s — — so, so). Since the flux J ~ zo, 
this condition guarantees zero flux at the drum wall. Equation for S is similar to Eq. ([l7|), but has an extra term due 
to rotation (compare 

d t 5 = n + d 2 s J (52) 

where J = ^ f(z ), compare Eq. (fl6|), and Q w const is proportional to rotation speed ui. The increase of the angle 
due to rotation is compensated by the flux of particle downhill described by the last term in Eq. (|5^) . 

We studied Eqs. (|5lf) , ([52|) numerically. The deep layer approximation is not valid near the edges of free surface, i.e. 
for s ±so. It results in anomalous growth of the angle 8 for s — > ±so. In order to prevent this spurious behavior 
we add regularization term — C(s)S to Eq. (|5^). The function £(s) was chosen as follows: ( — tanh(£o(so — l s l))i i- e - 
( — ► near the edges and C = 1 otherwise. In our numerical simulations we used Co = 0.2. We checked that the bulk 
behavior was not sensitive to the specific choice of the function ((s). 

Some of the results are presented in Figs. [I4|-|l6|. As seen in Fig. [l4|, for low rotation rate granular flow has a form 
of a sequence of avalanches separated by almost quiescent states (zq — > 0}. Surprisingly, the time behavior, especially 
in large drums, at low rotation rates fl is not strictly periodic, see Fig. H5l although one can distinguish well-defined 
characteristic time between the avalanches, like in Ref. |32|,|3|]. We think that this stochasticity in the form and the 
duration of avalanche events is related to the noise amplification (i.e. numerical noise). Since the avalanches are 
separated by long quiescent periods when zq is exceedingly small (zo could be as small as 10 -20 for £1 1), the slope 
of free surface 5 may "overshoot" the instability limit 5=1, and the system becomes susceptible to small fluctuations. 
These fluctuations trigger avalanches at random positions of the drum. 

For higher rotation speed we observed hysteretic transition to steady flow. In the steady flow regime (dtZo = dtS = 0) 
one finds from Eq. (m) 

J = (4 - s 3) (53) 

Using that J ~ Zq (see Eq. (fl6])), one immediately finds the dependence of the depth of fluidized layer zq on the 
position along the drum surface s: zq ~ f2 1//3 (s§ — s 2 ) 1 / 3 (this expression valid far from the edges, i.e |s| < so). The 
dependence of zq vs s for all values of s is shown in Fig. (ttq). As one sees from the figure, the dependence of zq vs s is 
practically symmetric with respect to the center of the drum, and Zq increases with the rotation rate Q, in agreement 
with experiments, see, e.g. j53|. The form of zq vs s appears to be consistent with recent experimental observations 
on flows in rotating drums Ref. [flij . 
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FIG. 14. Bifurcation diagram for rotating drum obtained from solution of Eqs. ( |5]] , ^ ) for fi — 0.2, (3 — 3.15, —so < s < so, 
so = 100. The symbols show zq at the center of the drum (s = 0) at the moments of time corresponding to dzo/dt = 0, • 
correspond to increase of fi, □ to decrease of fi. The arrows illustrate the hysteretic transition between stationary and avalanche 
flow. 
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FIG. 15. The width of fluidized layer zo vs time at the center of the drum in the regime of avalanche flow for three different 
values of fi, other parameters the same as for Fig. [l4|. 
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FIG. 16. Width of fluidized layer z vs s in the stationary flow regime for Q = 0.004; 0.0045; 0.005. 



V. SHEAR GRANULAR FLOWS AND GRANULAR FRICTION 



In this section we consider one-dimensional shear flow of granular matter placed between two parallel plates, one of 
which is moving with the velocity Vq (Couette flow). This flow (or rather Taylor-Couette flow between cylinders) has 
been studied in a number of recent experiments |^2|-|25| It was found that at small pulling speeds, the granular flows 
exhibit non-stationary stick-slip motion [Q. At higher pulling speeds, the flow becomes stationary. The velocity 
profile is typically exponential ]22| , |25| -p7[| in 2D experiments, but Gaussian in 3D Taylor-Couette geometry p3[ . We 
will show that these findings can be readily explained within our theoretical description. 



A. 2D problem 



In this section we neglect the effect of (bottom) immobile boundary and restrict our analysis by the case of planar 
shear flow in a semi- infinite layer of granular matter driven by a moving plate (see Fig.p/fl). The general model is 
reduced to Eq.(|l4]) combined with the constitutive relation 

°xz = ti— + pol z . (54) 
where v(z) is the horizontal velocity of the granular flow. 



FIG. 17. Schematic representation of a 2D granular flow experiment, 
is pulled via a spring with constant velocity Vq. 



Granular material is driven by a heavy top plate which 
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The relation for the control parameter S reads as 



(<7 x ,A&) 2 -# 



(55) 



where 4>i^ are tangents of the static and dynamic repose angles. Unlike the chute flow, the normal stress a zz is 
constant (=1), and the static component of the shear stress a xz is independent of the depth z. Here we neglect the 
weight of the sand itself since the weight of the top plate provides a much larger normal stress. The independence 
of the static stress a xz on z is in fact approximation. We assume that the "stress propagation time" , which is of the 
order of collision time To ~ \fdj 1 g is much smaller then any time scale in our problem, which is true for not too high 
shearing rates. 

The balance of forces requires a xz — a® z , which together with the constitutive relation ( |54| ) yields the expression 
for the shear velocity (one need to set ctq = o~ xz in order to satisfy the boundary condition at z — ► oo): 

v{z) = <r° xz I (1 - P )dz'. (56) 



These equations have to be augmented by the boundary conditions at z — > oo and z = 0. At z — > oo we require 
p — ► 1, i.e. the granular material is static, and at z = 0, we require the no-flux condition p z (0) = 0. In addition to 
that, we need a relation between the shear flow velocity and the shear stress near the surface. We propose that as 
the plate moves, the shear stress at the boundary is proportional to the difference between the displacement of the 
plate xq = Vot and the displacement of the grains immediately below the surface x (effective "Hooke's law"), 

°lz = 70o - x). (57) 

where 7 = const is proportional to the spring stiffness. Differentiating this equation with respect to time t yields the 
boundary condition 

&° xz = 7 (Vb - HO)), (58) 

Introducing scaled variables S = er° /(</>! - (j)\) 1/2 (J zz , So = <f>i/(<j>l - <j>\) 1/2 , V — 777, W — V r}/(4% - 
fil) 1 a zz' V( z ) = w ( z ) 7 7/('/ > 2 — (flY^^zzi we obtain the following set of equations, 

p = V 2 P + p(l-p)(p-S 2 + S 2 ), (59) 

V(z) =S f (1 - p)dz' (60) 

J — OO 

S = T(W - V(0)) (61) 

It is interesting to trace the connection of our boundary condition (^Tj) to the Maxwell stress relaxation condition 
for visco-elastic fluid p6| 

d °xz 1 ^ du xz 

it + = E ^r (62) 

where E — const is the "shear modulus" and u xz is the strain tensor, and r is some characteristic relaxation time. It 
can be expected that r is a function of order parameter and diverges in solid state, p — ► 1, so if we take 

r = r/(i- P ), (63) 

and integrate Eq. ( |62] ) over z from —00 to 0, we get Eq.(|6l|). 

We integrated the equations (|59|),(60),(|6"f|) numerically using finite difference method. The main control parameter 



of this model is the normalized velocity W. At large W we obtained a stationary near-surface shear flow with a profile 
shown in Fig. [l^. In fact, this stationary distribution of the order parameter p coincides with the exact solution Eq. 
(fH|), in which S should be replaced by S 2 — S 2 . 
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At small velocities W — > the model exhibits relaxation oscillations, reminiscent of the normal dry friction between 
two solids. The stress a xz grows almost linearly with no flow until it reaches a certain threshold value after which the 
near-surface layer fluidizes, and the ensuing shear flow relieves the accumulated stress. After that the layer "freezes" 
again, and the process repeats (see Fig. |19|) . 
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FIG. 19. Relaxation oscillations of the shear stress o xz (a) and the near-surface velocity V(0) (b) for 
So = 0.1, T = 0.01, W = 2. 

Fig. [2(] depicts the bifurcation diagram illustrating the transition from the stationary shear flow at large W to the 
regime of relaxation oscillations at small W . As can be seen, the transition is subcritical with hysteresis (similar to 
that occurred for rotating drum), as the oscillations always occur with the finite amplitude. Similar abrupt transition 
from oscillations to steady sliding was found in experiment Ref. |^4j . 

Furthermore, we explored the dependence of the shear stress S on the pulling speed W . For the stationary flow 
regime it can be done fully analytically using the exact solution Eq. (21 ) for the order parameter p. Simple expression 
can be obtained at the large velocity limit (W 3> 1) corresponding to 5 — * 1/2. In this case one derives from Eqs. 
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w = s 



(1 - p)dz' 




(64) 



As one sees from Eq. (|6J) , with increase of the pulling velocity W the shear stress S monotonically decreases and 
approaches the value S = \/l/2 + Sq, in agreement with experimental results of Ref. |25f| where it was also found 
that the shear stress slightly decreases with the increase of shear rate and approaches some equilibrium value. 

The numerically obtained dependence of S vs W for arbitrary values of W is shown in Figure As seen from the 
Figure, indeed there is only a weak dependence of S on W, up to rather small value of W sw 0.5. 




FIG. 20. Bifurcation diagram for the transition from stationary shear flow to relaxation oscillations. The dots in this plot 
depict the extrema of V(t) as a function of the pulling speed W for So = 0.1, F = 0.01. Periodic oscillations coexist with steady 
sliding for 2.6 < W < 3.6. 




FIG. 21. Dependence of the normalized shear stress S on the pulling speed W for So = 0.1, T — 0.01 (circles). Non-uniqueness 
of this function is a result of the hysteretic transition from stick-slip to continuous sliding. Solid line corresponds to solution 
( pi] ) for stationary shear flow regime. 

In recent paper Ej ], a model for the granular friction has been proposed which is based on similar ideas of a 
phase transition in the granular medium underneath of a moving surface. Within this model, oscillations in the form 
of stick-slip motion can be described, however the model does not described the observed transition from stick-slip 
motion to the steady motion with increasing of the pulling velocity. The significant difference between the model 
jlTj and our model is that the former does not address the spatial inhomogeneity of the fluidized layer, and thus the 
order parameter dynamics is described by an ordinary differential equation. Second, the control parameter in the 
order parameter equation |47| (see also ]4q] ) is a function of the sliding velocity and not the applied stress, which 
in our opinion is not physical. Indeed, the transition to a fluidized state should be determined by a yield condition 
which is naturally defined via components of the stress tensor. Although in the dynamic friction problem the sliding 
velocity and the shear stress are related, we believe that the motion of the granular material is the consequence of 
the fluidization transition rather than the reason for it. 

An alternative approach leading to the exponential decay of velocity in shear flows was developed in Refs. p5|. 
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Using traditional hydrodynamic equations coupled to the equation for the granular temperature, the authors repro- 
duced experimentally behavior. However, in order to explain the exponentially small velocity tail far away from the 
moving plate, highly-nonlinear viscosity with singular dependence on the density was introduced ad hoc. That model 
successfully describes the shear flow driven by a moving plate, however, it fails to describe the transition between 
solid and fiuidized states which is the hallmark of the granular dynamics. Our model, on the contrary, is applicable 
to description of both flowing and static regime. 



B. 3D problem 



In this Section we consider the 3D shear flow structure between two vertical plates one of which is moving with 
respect to the other (see Fig. ^||). This geometry is inspired by the recent experiment by Mueth et al. |2i| in which 
the structure of the granular shear flow was studied in a long vertical Taylor-Couette cell. They found a significant 
deviation in the shear flow profile from the simple exponential profile observed in earlier 2D experiments |£2 24 2j| 
and successfully reproduced by our theory (see the previous Section). 




FIG. 22. Schematic representation of a 3D shear experiment. A slab of granular material is sandwiched between two vertical 
plates, one is moving with the speed Vo at y — and one immobile at y — — 2L 



Ref. |23| gives strong evidence for the Gaussian (v ~ exp(— const x (r — r^) 2 ) behavior of the velocity profile near 
the outer wall. As we show below, this feature can be attributed to essentially three-dimensional geometry of the 



experiment in contrast to that of Refs. p2|j2 
in sufficiently long vertical cylinders filled wit 



It follows directly from the fact that the normal stress (pressure) 
h dry granular materials, saturates and does not depend on height of 
the cylinder (as in the celebrated Jansen picture of silo, see ||). 

Indeed, consider the distribution of stresses in infinite layer of grains in contact with 2 vertical walls (gravity is 
directed along the z-axis) at the points y = and y — —2L, see Fig. ^2|. From the projection of force on z axis we 
obtain the condition 



'yz,y 



-g- 



(65) 



Since we assume that the diagonal component a zz does not depend on depth, i.e. z-coordinate, and there is no 
cc-dependence of stresses, we obtain that the weight is supported by the tangential stress a yz = —g(y + L). Due to 
shearing we will also have tangential stress o~ xy — const, as in the two-dimensional case (with notational difference 
that y now indicates the direction normal to the walls). 

Now, we need to relate the tangential stresses to the control parameter S using the Mohr- Coulomb condition in 
three dimensions. It is well known that the tangential and normal stresses a T and o~ n lie within the shaded area 
limited by the Mohr's circles built upon the major, intermediate, and minor principal stress values <ti_3. (see Fig.p3|). 
As it is clear from Fig.^, the maximum value of the ratio <j) = o~ T / 1 o~ n occurs at the tangential point A at which 



o~i — 0~3 



(66) 



The major and minor principal stresses are determined as eigenvalues of the shear stress tensor o~ij. The control 
parameter 8 is related to (j> via Eq.(||), so in this case 



>i -<r 3 )2 
4ctict 3 



(67) 
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FIG. 23. Mohr's circle for the three-dimensional shear stress tensor M. 

For simplicity we assume that all diagonal components of the stress tensor are equal (like in fluid): a xx = a yy — 
a zz = P = const. This assumption makes the calculation much simpler, although qualitatively similar behavior 
is expected in general case. Here P has the meaning of pressure on the bottom of the cell and from dimension 
consideration one concludes that P cogL, where Co > 1 is a constant which depends on surface wall roughness etc. 
"291. One obtains the eigenvalues 



g 2 (y + L) 2 + a 2 



0-1 = P- 

<J 2 =P 



(68) 



Thus, one derives 



S(y) = 



g 2 (y + L) 2 + a : 



xy 



p2-g*(y + L) 2 -a 2 



(69) 



Here, as before, a xy characterizes shear rate and is proportional to Vq. Therefore, the control parameter 5 in 3D has 
an explicit y-dependence. In contrast, in 2D geometry we had simply S = const. 

Let us now evaluate the decay rate of the order parameter. Linearizing Eq. (|0|) near p = 1 we obtain 



Pyy ~ Pi 1 - S (V)) = 



(70) 



We focus on the solution near the wall y — 0. If \y/L\ <C 1, we can apply WKB approximation and seek the solution 
in the from p ~ exp($(y)) into Eq. ((?o|). Then we derive 



$2 



1 - %) 

Therefore, for < \y \ <C L and also for a xy /P <C 1 we derive 



exp 



y/1 - S(y>)dy' 



exp 



V^-S(0)y 



S'(0) 



Vi - s(o) 



V 2 + 0(y 3 ) 



dy' 



(71) 



(72) 



Thus, one sees that in 3D p(y) possesses Gaussian correction term which is absent in 2D because 6' = in the 2D 
shear flow. Let us point out that, to the first order, the coefficient in front of y 2 does not depend on the shear rate 
in agreement with experiment (23|. 



VI. CONCLUSIONS 



We developed a continuum theory of partially fluidized granular flows. This theory is based on a combination of 
the mass and momentum conservation laws with an equation for the order parameter describing the transition from 
the static to flowi ng reg ime. In this sense, our theory goes beyond the hydrodynamical description of dense granular 
flows, see e.g. p9|-pl 24 . The order parameter which is a crucial variable in our theory, can be interpreted as a portion 
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of the static contacts among particles in a small volume within the granular system. This characteristic is difficult 
to measure in physical experiments, but can be extracted from molecular dynamics simulations. Phenomenological 
parameters in our model can be obtained from comparison with molecular dynamics simulations and experiments. In a 
certain limit our model can be reduced to two coupled equations for the depth of fluidized layer and local angle, which 
resemble BCRE model, however differ from it in detail. In particular, our model has intrinsic hysteretic behavior 
absent in BCRE model. 

It may appear surprising that our theory has the characteristic length of the order of the grain size d. Usually, the 
hydrodynamic description is valid on the scale much larger then d. We believe that in case of dense granular flows, the 
continuum (but non-hydrodynamic) description is possible because in the flow regime the granular system "samples" 
all possible states in the configuration space, thus providing "self-averaging" . Moreover, our model exhibits "critical 
slowdown" for 8 — ► 1. In this case the decay length (or, pursuing the analogy with equilibrium critical phenomena, 
the "coherence length") d s — d/y/1 — 8 diverges at the critical point. Thus, rigorous derivation of order parameter 
equation can be anticipated in the vicinity of point of spontaneous fluidization 8 = 1. 

Our order-parameter model captures many important aspects of the phenomenology of chute flows observed in recent 
experiments [|l^- pl| , ^6| , p7| , including the structure of the stability diagram, triangular shape of downhill avalanches at 
small inclination angles and balloon shape of uphill avalanches for larger angles. It provides an adequate description 
of granular flows in a 2D rotating drum and in Couette geometry. In particular, we found the experimentally observed 
features such as periodic oscillations of the shear stress and flow velocity at low rotation rates and transition to steady 
flow at higher rates. For shear cell experiment our model gives rise to the exponential velocity profile in two dimension 
and Gaussian correction to the profile in three dimensions. 

We believe that our model can be applicable to other granular flows and can be generalized to binary mixtures of 
granular materials, wet and cohesive granular materials and granular flows with additional vibration. Another chal- 
lenging project is to derive the order parameter description from some sort of "microscopic" theory of granular flow, 
in analogy to the theory of superconductivity, where the order parameter equation was first proposed phenomenolog- 
ically by Ginzburg and Landau |54| , and later derived from the microscopic theory of superconductivity by Gorkov 
and Eliashberg |5q| . 

We thank Dan Howell, Saturo Nasuno, Doug Durian, Pierre-Gilles de Gennes, Jerry Gollub, Leo Kadanoff, Bob 
Bchringer and Adrian Daerr for useful discussions. This research is supported by the Office of the Basic Energy Sciences 
at the US Department of Energy, grants W-31-109-ENG-38, DE-FG03-95ER14516, and DE-FG03-96ER14592. 
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APPENDIX A: THE TRANSITION LINE BETWEEN TRIANGULAR AND UPHILL AVALANCHES: 

INFINITE VISCOSITY LIMIT 



The transition line between triangular and uphill avalanches can be found analytically in the limit of infinite viscosity. 
In this case Eq. ( p"l| ) yields trivial solution h = ho = const and Eq. ([l4j) becomes independent. Within a single order 
parameter equation, the velocity gap disappears, and the transition line corresponds to V = 0. 
We represent the solution for the order parameter in the from 



p(x, z, t) — p(x + Vt, z) 
Substituting the ansatz (Al) in Eq. |l4| one obtains 

V P , = V 2 p-p(l-p)(S-p) 



For V = Eq. (A.2) is reduced to 



V 2 po-po(l-po)(<5-po) =0 



(Al) 
(A2) 
(A3) 



The solution to Eq. ( |A3j ) exists only for some specific value of S = So for each ho. This value can be obtained from the 
solvability condition. Since d x po is the solution of the linearized problem, the solvability condition is obtained with respect 
to this solution. Multiplying Eq. (A3) by d x po and performing integration over x and z, one obtains 



2G 



dx 



dzd x p (V 2 po - p (l - po)(5o - po)) = 



Using integration by parts, we derive 



dx j dz-^- ((d z p ) 2 + ipo - |(*o + l)po + <5po) = 



(A4) 



(A5) 



which leads to 



dz ({dzpo) 2 + ipo - |(<5o + l)po + Spa 



h (25-l) 
6 



= 



For x — ► oo the solution po converges to a pure one-dimensional solution with the first integral 

- (d z po) 2 + ipo ~ g(<5o + l)po +$po = C = const 



where C : 



2P0 



j (5o + l)po + <5po at z = 0. Therefore, the expression (A6) can be brought to the form 



J dz Qpo ~ |(^o + l)po + <5po) 



gh _fto(M-l) =0 



(A6) 



(A7) 



(A8) 



Setting h — ho and solving Eq. (AS) for each /i, one find the dependence 80 vs ho- This dependence is shown in Fig. [ll], 
dotted line. The infinite viscosity limit gives the lower bound of uphill avalanches, however, as one sees in the Figure, this 
limit is rather close to the experimental data. 

For /1 > 1 one can derive an estimate for h vs 5. In this limit, the solution is given by the front solution Eq. (|34|). Since 
for 8 > 1/2 the fluidized state invades the solid state, the front travels toward t he b ottom and stops at the distance 
Az ~ log(<5 — 1/2) ~ O(l) from the bottom. Thus, for /i> 1 one obtains from Eq. (A8) (taking into account C — > 0) 



h • 



5-1/2 ' 



(A9) 



Let us now consider the case h ~ 0(1). In this case the transition line h(5) can be obtained analytically from the single- 
mode approximation Eq. (pi]). To demonstrate that, we first find the position of the line V — in the (5, h) plane for 
ol — > 0. Eq. ( p4| ) has a stationary front solution connecting two outer fixed points A\ and A3 of Eq.(^iJ) if its free energy is 
symmetric, i.e. the roots j4i,2,3 of equation \(h)A+ 3^ 
It gives rise to expression 



8(2 S) A 2 - fA 3 = satisfy the symmetry condition Ai = 0, A s = 2A 2 . 



7T 2/16 

4/? + 3 



0. 



From Eq. (A10) we obtain 



1 + 



.(2-5)) ; 



(A10) 



(All) 



The dependence h(S) agrees with the corresponding dependence obtained from the analysis of full Eq. (|l^ ) (see Eq. JA8|) ) 
within the line thickness up to <5 ~ 0.6. Also, along this line there is an exact expression for the front solution 

(A12) 



A = A x (l + tanh(-s/3/8Aoox)j 



where Aoa is given by (B 



APPENDIX B: UPHILL FRONT VELOCITY IN THE LARGE VISCOSITY LIMIT 

For finite a we look for the solution in the form 



A = Aq[x) + eAi 
h = ho + ehi (x) 
V = eVi 

S = 5* + eS\ — f3ed x hi 



(Bl) 
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where e = y/a. Substituting the ansatz (Bl) into Eqs. (p4|), (pq), we obtain in the first order in e: 

8 * - 2 - 



LAi 



Vihi 



ttAi - /Jfc! (— A 2 - A ) 



-h A 



(B2) 
(B3) 



where L is linearized Eq. ( J24[ ) in the vicinity of An at h = ho,5 = 5* . 

Eq. (B2) has bounded solution if the r.h.s. of Eq.(B2) with hi expressed from Eq.(B3) is orthogonal to the zero eigenmode 
Ai = d x Aq. This solvability condition yields 



Viai + — 5i«3 = 
Vi 



where 



(B4) 



ffll = 



a 2 = 



a 3 



(^Ao) 2 ^ = v/273AL 



/3^(a.A ) 2 (^A 2 , - A ) + 7T 2 /2A f 2 AA 



dx = A; 



4tt 2 16/3 V^gA 2 ^ _ /3^6fegA c 
3 + 15tt 3 



(A - ^-Ag)^A = 2AL 

37T 



^A 3 



(B5) 



Returning to the original notation, we obtain from Eq. (B4) 



V-Sd 1 V + ad 2 = 



(B6) 



where di = 0,3/01, d2 = 0,2/0,1 and 5 = 5 — 5*. 
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